#include <kvs/TetrahedralCell>
#include <kvs/UnstructuredVolumeObject>
#include <kvs/UnstructuredVolumeImporter>
#include <kvs/StructuredVolumeObject>
#include <kvs/ValueArray>
#include <kvs/Type>
#include <kvs/glew/RayCastingRenderer>
#include <kvs/glut/Screen>
#include <kvs/glut/Application>


int nx = 480;
int ny = 720;
int nz = 120;
double valid_opacity = 6.0/256.0;

void CalcBoundingBox( kvs::TetrahedralCell<float>* cell, kvs::Vector3f& min_c, kvs::Vector3f& max_c )
{
    kvs::Vector3f c[4];
    for ( int i = 0; i < 4; ++i )
    {
        c[i] = cell->vertices()[i];
    }

    min_c = c[0];
    max_c = c[0];
    for ( int i = 1; i < 4; ++i )
    {
        for ( int j = 0; j < 3; ++j )
        {
            min_c[j] = kvs::Math::Min( min_c[j], c[i][j] );
            max_c[j] = kvs::Math::Max( max_c[j], c[i][j] );
        }
    }
}

kvs::StructuredVolumeObject* FloatImporter( std::string filename, size_t nx, size_t ny, size_t nz )
{
    FILE *rawfile = fopen( filename.c_str(), "rb" );
    size_t length = nx * ny * nz;    
    
    std::cout << length << std::endl;
    
    kvs::AnyValueArray values;
    float* buf = static_cast<float*>( values.allocate<float>( length ) );
    fread( buf, sizeof( float ), length, rawfile );
    fclose(rawfile) ;
    
    kvs::Vector3ui resolution( nx, ny, nz );
    kvs::StructuredVolumeObject* object = new kvs::StructuredVolumeObject(
                                                                          resolution,
                                                                          1,
                                                                          values
                                                                          );
    object->updateMinMaxCoords();
    object->updateMinMaxValues();
    
    std::cout << "Loading of jet file succeed" << std::endl;
    std::cout << *object << std::endl;
    
    return( object );
}

int main( int argc, char** argv )
{
    //kvs::glut::Application app( argc, argv );
    kvs::StructuredVolumeObject* ori_volume = FloatImporter( argv[1], nx, ny, nz );;
    kvs::UnstructuredVolumeObject* eva_volume = new kvs::UnstructuredVolumeImporter( argv[2] );
    float max_value = (float)ori_volume->maxValue();
    std::cout << max_value << std::endl;
    
    float* ori_values = new float[nx * ny * nz];
    ori_values = (float*)ori_volume->values().pointer();
    
    size_t NCells = eva_volume->ncells();
    kvs::TetrahedralCell<float>* cell = new kvs::TetrahedralCell<float>( eva_volume );
    
    kvs::ValueArray<float> eva_values;
    eva_values.allocate( nx * ny * nz );
    
    for ( size_t count = 0; count < NCells ; count++ )
    {
        cell->bindCell( count );
    
        kvs::Vector3f minf, maxf;
        kvs::Vector3i mini, maxi;
        CalcBoundingBox( cell, minf, maxf );
        
        for ( int a = 0; a < 3; ++a )
        {
            mini[a] = kvs::Math::Floor( minf[a] ); 
            maxi[a] = kvs::Math::Ceil( maxf[a] ); 
        }
        
        if ( mini.x() < 0 || maxi.x() > nx ) std::cout << "x out" << std::endl;
        if ( mini.y() < 0 || maxi.y() > ny ) std::cout << "y out" << std::endl;
        if ( mini.z() < 0 || maxi.z() > nz ) std::cout << "z out" << std::endl;
        
        for ( int k = mini.z(); k <= maxi.z(); k++ )
            for ( int j = mini.y(); j <= maxi.y(); j ++ )
                for ( int i = mini.x(); i <= maxi.x(); i ++ )
                {
                    cell->setGlobalPoint( kvs::Vector3f( i, j, k ) );
                    kvs::Vector3f local = cell->localPoint();
                    if ( !(local.x() < 0 || local.y() < 0 || local.z() < 0 || 1 - local.x() + local.y() + local.z() < 0 ) )
                    {
                        int index = i + j * nx + k * nx * ny;
                        eva_values[ index ] = cell->scalar();
                    }
                }
    }
    
    //Calculate the error
    float error;
    float average_error;
    size_t count = 0;
    size_t count_error = 0;
    
    double valid_value = max_value * valid_opacity;
    std::cout << "valid value: " << valid_value << std::endl;
    for ( long long i = 0; i < nx * ny * nz; i ++ )
    {
        if(ori_values[i] >= valid_value)
        {
            error += fabs( (eva_values[i] - ori_values[i]) / max_value ) ;
            count ++;
            if( eva_values[i] - ori_values[i] ) count_error++;
        }    
    }
//    average_error = error/(nx * ny * nz);
    average_error = error / count;
    
    std::cout << "count: " << count << std::endl;
    std::cout << "count_error: " << count_error << std::endl;
    std::cout << "error rate: " << (float)count_error/count << std::endl;
    
    std::cout << "whole error: " << error << std::endl;
    std::cout << "The average error is: " << average_error << std::endl;
    

}
